
template<
  int deal_num_per_block, int block_sz>
__global__ void reduce_sm75(
  const int m, const float* __restrict__ src,
  float* __restrict__ dst)
{
  __syncthreads();
  const int tx = threadIdx.x;
  const int bx = blockIdx.x;
  const int start = bx*deal_num_per_block;
  const int end = start+deal_num_per_block;
  __shared__ float stmp[block_sz];

  float sum = 0.f;

  for(int id = start+tx; id < end; id+=block_sz){
    if(id < m)
      sum += src[id];
  }
  stmp[tx] = sum;
  __syncthreads();

  for(int stride = block_sz/2; stride > 0; stride >>=1){
    if(tx < stride){
      stmp[tx] += stmp[tx+stride];
    }
    __syncthreads();
  }
  if(tx == 0){
    dst[bx] = stmp[0];
  }
}

int reduce(
  const int m, void *src, void *dst,
  cudaStream_t cuStream = nullptr)
{
  using SDT = float;
  using DDT = float;

  SDT *d_src = (SDT*)src;
  DDT *d_dst = (DDT*)dst;
  const int deal_num_per_block = 1024;
  const int block_sz = 256;

  dim3 g_sz(DIV_UP(m, deal_num_per_block)), b_sz(block_sz);
  reduce_sm75<deal_num_per_block, block_sz><<<g_sz, b_sz, 0, cuStream>>>(m, d_src, d_dst);
  return 0;
}